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Abstract 

A multiplication operator on a Hilbert space may be approximated with finite sections by choosing an 
orthonormal basis of the Hilbert space. Nonzero multiplication operators on L 2 spaces of functions are 
never compact and then such approximations cannot converge in the norm topology. Instead, we consider 
how well the spectra of the finite sections approximate the spectrum of the multiplication operator whose 
expression is simply given by the essential range of the symbol (i.e. the multiplier). We discuss the case of 
real orthogonal polynomial bases and the relations with the classical Fourier basis whose choice leads to 
well studied Toeplitz case. The use of circulant approximations leads to constructive algorithms working 
for the separable multivariate and matrix-valued cases as well. 

Keywords: Multiplication operator, orthogonal polynomials, Fourier basis, Toeplitz (and Generalized 
Locally Toeplitz) sequences, symbol. 


1 Introduction 

This note is in some sense a consequence of the intriguing MathSciNet Revue by Albrecht Bottcher of a paper 
by Morrison [13] and, of course, of the intriguing paper itself. Briefly, if 0 is a bounded function defined 
on a compact set I< of R d , d > 1, consider the multiplication operator M[cj>] : L^,(K) —» L^(K) defined as 
M[</>]( f) = <j>f, w suitable weight function. It is known that the spectrum is given by the essential range of <j>: 
now suppose that we have only a finite number of coefficients (Mjf [</>]) t . = {M[<f>]ej,ei), i, j = 0,..., N — 1, 
with {ey} denoting an orthonormal basis of the question is about the reconstruction of the multiplier <fi 
from the spectra of Mn[4>]. For reconstruction we mean the convergence of the finite sections spectra to the 
the essential range of the symbol <fi. More in general, we are interested in understanding as much as possible 
about (f> only using the entries of the matrices for large but finite N. 

Indeed the problem posed is a classical one (a beautiful historical account can be found in [13]). Here 
the idea is to discuss how the case of the choice of a general real orthogonal basis on I\ = [—1,1] can be 
reduced to the Fourier case and therefore to the Toeplitz case (see [4, 5, 11] and references therein for an 
encyclopedic coverage from three different angles) and how the latter can be reduced to the circulant case. 
Circulants (see [6]) are normal matrices since they can all be diagonalized by the same unitary transform. 
Further the transform is the celebrated discrete Fourier transform (DFT) for which a stable and extremely 
efficient algorithm exists (the Fast Fourier Transform i.e. FFT, see [33]). Therefore the general case can be 
translated into a problem of (asymptotic) structured numerical linear algebra for which an accurate solution 
can be determined with a low computational cost (here for low cost we mean 0{N\og{N)) arithmetic 
operations i.e. the asymptotic cost of a generic FFT). Moreover, the restriction on the boundedness of (f> 
can be suppressed and, more precisely, a related symbol <j> (more specifically 4>{x) = cf>(x)w(x)V 1 — x 2 ) has 
to be supposed just Lebesgue integrable: in this case, the operator M[cj)\ can be unbounded and has to be 
defined on a different domain. Multidimensional block generalizations (for multiplication operators having 
a matrix-valued multivariate function as multiplier) are also available thanks to the rich theory built in the 
finite dimensional case in recent years. The paper contains three more sections: Section 2 is devoted to 
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linear algebra tools; in Section 3 we discuss the solution to our problem and we give a brief account on 
separable multivariate and matrix-valued generalizations; Section 4 is concerned with open questions and 
final remarks. 


2 Notation from asymptotic linear algebra 

First we introduce some notations and definitions concerning general sequences of matrices. For any function 
F defined on C and for any matrix A n of size d n , with eigenvalues A j(A n ) and singular values <jj(A n ), 
j = 1,..., dn, by the symbols E (T (F, A n ) and T,\(F, A n ) we denote the means 

j-J^FiaMn)}, ±.f^F[\ j{ A n )], 

n j=1 " 3 = 1 

and by the symbol || • || the spectral norm i.e. ||X|| is the maximal singular value of the matrix X (see [1]). 
Furthermore || • || p indicates the Schatten p norms, p £ [l,oo) defined as 

\\A n \\ p p = Z a (\ • | p ,A n )-d n . 

The Schatten oo (p = oo) norm is exactly the spectral norm (for a unified treatment of these norms refer to 
the beautiful book by Bhatia [1]). Moreover, given a sequence {A n } of matrices of size d n with d n < d n+ 1 
and given a /z-measurable function g defined over a set K equipped with a a finite measure p, we say that 
{ A n } is distributed as (g, I\, p) in the sense of the singular values (in the sense of the eigenvalues) if for any 
continuous F with bounded support the following limit relation holds 

(1) lim T, a {F,A n ) = —[ F(\g\)dp, ( lim T, x (F,A n ) = — l — [ F(g)dp\ . 

n^oo p(K)J k \n ^oo p(K) J K J 

In this case we write in short {A n } ~ <T ( g , K, p) ({ A n } ~a {g, K, p)). An interesting connection between the 
notion of distribution and the Schatten p norms is given in the following Lemma. 

Lemma 2.1 Assume that {A„} ( g,K,p) and that ||B ra || p = o(dl/ p ), A ni B n both of size d n , and p £ 

[l,oo]. Then it holds 

(2) {B n } (0, K, p) and {A n + B n } ( g , K, p). 

Moreover, if all the involved sequences are Hermitian and {A n } ~a (g,K,p), then (2) holds true with ~ CT 
replaced by 

Proof. The tools for the proof in the case of p = 2 can be found in [31]. Here we treat the general 
case by using analogous ideas. For p = oo and ||f? n || = o(l) the proof is trivial by standard perturbation 
arguments (see e.g. [1, 34]). Therefore we focus our attention on the case where p £ [l,oo). Indeed, from 
the assumptions on {B„j with p £ [l,oo), for every e > 0, we have 

d n 

C(n) = ||£n||£ = $>?(£„) 

3 = 1 

crj(B n )>e 

* E ^ 

aj(B n )>e 

= e p #{a j (B n ) > e} 

with C(n) = o(d n ). Therefore the cardinality of the singular values bigger than e is bounded from above by 
C(n)/e p = o(d n ). Since e > 0 is arbitrary, by direct check, it follows that { B n } (0, I\ , p). Furthermore, by 

exploiting the singular values decomposition of B n , we can write B n as L n (e) and R n (e) where ||L ra (e)|| <x) < e 
and the rank(l? n (e) < C(n)/e p = o(d n ). More precisely, in the previous lines we have proved that the 
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cardinality of the singular values of B n bigger than e is bounded from above by C{n)/e p = o{d n ). Now 
from the SVD decomposition (see e.g. [1]) there exist U n and V n unitary matrices and D n diagonal matrix 
(containing the singular values of B„ sorted non decreasingly) such that 

B n = U n D n V n . 

At this moment take D n (>) the matrix containing all the entries bigger than e of D n (in the same position 
as D n ) and D n (<) the matrix containing all the entries at most equal to e of D n (in the same position as 
D n ). Therefore D n = D n {>) + D n (<) with 

||-Dn(<)|| < e, rank(D„(>)) < C(n)/e p = o(d n ). 

Finally since U n and V n are unitary we have 

\\U n D n (<)V n \\ = ||D n (<)|| < e, rank(U n D n (>)V n ) = rank(£>„(>)) < C(n)/e p = o{d n ), 

and B n = U n D n (<)V n + U n D n (>)V n . The statement is proven by putting L n (e) = U n D n (<)V n and 
R n (e) = U n D n (>)V n 

Consequently, by using e.g. Proposition 2.3 and Remark 2.1 in [19], from the hypothesis {A n } (g, K, g) 
we deduce {A n + B n } (g, I\, g). The case of the eigenvalues for Hermitian matrices A n and B n is identical 
and it is not repeated here. • 


2.1 How to use spectral distribution 

We show how the notion of distribution can be used for the reconstruction of the symbol when the eigenvalues 
(or singular values) are known. More precisely, the subsequent Theorem 2.1 demonstrates that { A n } 

(, g,K,g ) (or {A n } ( g,K,g )) and the knowledge of the eigenvalues of {A n } (or singular values of { A n }) 

imply that many facts on the symbol g can be constructively recovered. 


Definition 2.1 Given the g measurable function g defined on K with g being a a finite measure supported 
on K, the (essential) range of g is given by the points s € C such that, for every e > 0, the measure of the 
set {t £ D : g(t) £ D(s,e)} is positive with D(s,e) = {z £ C : |z — s| < e}. The function g is (essentially) 
bounded if its essential range is bounded. Finally, if g is real-valued then the (essential) supremum is defined 
as the supremum of its range and the (essential) infimum is defined as the infimum of its range. 


Definition 2.2 A sequence {A n } (A n of size d n ) is properly (or strongly) clustered at s £ C in the eigenvalue 
sense, if for any e > 0 the number of the eigenvalues of A n not belonging to D(s,e ) = {z £ C : |z — s| < e} 
can be bounded by a pure constant q e possibly depending on e but not on n. Of course if every A n has, at least 
definitely, only real eigenvalues, then s has to be real and the disk D(s, e) reduces to the interval (s — e, s + e). 
Furthermore, a sequence {A n } (A n of size d n ) is properly (or strongly) clustered at the nonempty closed 
set S C C in the eigenvalue sense if for any e > 0 the number of the eigenvalues of A n not belonging to 
D(S, e) = u D(s,e) can be bounded by a pure constant q e possibly depending on e but not on n and if every 


sGS 


A n has, at least definitely, only real eigenvalues, then S has to be a nonempty closed subset of R. The term 
“properly (or strongly)” is replaced by “weakly ” if q e is a possibly unbounded function of n with q e (n) = o(d n ) 

Qe f 77<) 

(i.e. lim —— = 0). Finally, the above notions are in the singular value sense if the term “eigenvalue” is 


replaced by “singidar value”: of course s has to be a real nonnegative number and S has to be a subset of 
nonnegative numbers. 


Definition 2.3 A sequence {A n } (A n of size d n and with spectrum T, n ) is strongly attracted by s £ C if 


lim dist(s, E ra ) = 0 

n —>oo 

where dist(X, Y) is the usual Euclidean distance between two subsets X and Y of the complex plane. Fur¬ 
thermore, let us order the eigenvalues according to its distance from s i.e. |Ai(A n ) — s| < |A 2 (A n ) — s| < 
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■ ■ ■ < |A d„(A n ) — s\. We say that the attraction is of order r(s) £ N, r(s) > 1, fixed number independent of 
n, if 

lim |A r(s )(A„) - s\ = 0, liminf |A r ( s ) + i(A n ) - s| > 0. 

n—> oo n—>oc 

The attraction is of order r(s) = oo if 

lim |A j(A n ) - s\ = 0 

n—> oo 

for every fixed j independent of n. Finally, the term “strong or strongly” is replaced by “weak or weakly” if 
every symbol lim is replaced by lim inf. Finally, the above notions are in the singular value sense if the term 
“eigenvalue ” is replaced by “singular value”, T, n is replaced by the set of the singular values, and, of course, 
the value s is a real nonnegative number. 

Notice that writing {A n } ( g , K , p) with g constant function equal to s £ C is equivalent to write 
that {A n } is weakly clustered at s in the eigenvalue sense. Analogously, writing {A n } (<?, K , p) with g 
constant function equal to s £ R, s > 0, is equivalent to write that {A n } is weakly clustered at s in the 
singular value sense. 

The notions previously introduced are intimately related as emphasized in the subsequent theorem which 
is explicitly given only for the eigenvalues (the singular value version is obvious and is shortly sketched). 

Theorem 2.1 Let {A n } be a matrix sequence with A n having size d n and let g be a p-measurable function 
defined on K with p being a finite measure supported on I\. Consider the following statements: 

a ) { A n} (. g,K,p ); 

b) the (essential) range of g is a weak cluster for {A„} in the eigenvalue sense; 

c) the (essential) range of g strongly attracts the eigenvalues of {A n }; 

d) any point s of the (essential) range of g strongly attracts the eigenvalues of {A n } with order r(s ) = oo. 

e) given s £ C, e > 0, if the cardinality of the eigenvalues of A n belonging to D(s , e) divided by d n tends to 

a positive value, then s belongs to the (essential) range of g within an error of at most e; 

f) given s £ C, e > 0, if the cardinality of the eigenvalues of A n belonging to D(s,e) divided by d„ tends to 

a zero, then s cannot belong to the (essential) range of g. 

Then a) implies b), c), d), e), and f). Finally, the above implications hold in the singular value sense if 
the term “eigenvalue” is replaced by “singular value”, g is replaced by \g\, and of course the value s is a real 
nonnegative number. 

Proof. The first three implications are proven in Theorem 2.7 of [9]. For the other two see e.g. Section 4 
in [17]. • 

In the rest of the paper, with regard to relationships (1), the symbol p is suppressed for the cases under 
study (Toeplitz sequences, Generalized Locally Toeplitz sequences, Circulants etc.) since the measure will 
always coincide with the standard Lebesgue measure on R d for some positive integer d. 

2.2 Toeplitz matrix sequences 

Let rn{-} be the Lebesgue measure on R d for some d and let / be a d variate complex-valued (Lebesgue) 
integrable function, defined over the hypercube Q d , with Q = (—7r, n) and d > 1. From the Fourier coefficients 

of / 

(3) fj = ^ ndx f /(s)exp(-i(j,s)) ds, i 2 = -1, j = (ji, • • -, jd) £ Z d 

J JQ d 

with (j,s) = n = {n\,... ,Ud) and N{n) = n\---nd, we can build the sequence of Toeplitz 

matrices {T n (/)|, where T n (f) ={/j-i}™ J=1 T € C N ( n )x N ( n ), 1 T = (1,..., 1) g N d is said to be the Toeplitz 
matrix of order n generated by /. Furthermore, throughout the paper when we write n —> oo with n = 
(ni,..., Ud) being a multi-index, we mean that mini<j<d nj —> oo. 
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The asymptotic distribution of eigen and singular values of a sequence of Toeplitz matrices has been 
thoroughly studied in the last century (for example see [4] and the references reported therein). Here we 
report a famous Theorem of Szego [10], which we state in the Tyrtyshnikov and Zamarashkin version [32]: 

Theorem 2.2 If f is integrable over Q d , and if {T n (/)} is the sequence of Toeplitz matrices generated by 
f, then it holds 

(4) {T n (f)}~Af,Q d )- 

Moreover, if f is also real-valued, then each matrix T n (f) is Hermitian and 

(5) {T n (f)}~x(f,Q d ). 

This result has been generalized to the case where / is matrix-valued (see, for example, [28, 16] and 
Subsection 3.3) so that the matrices T n (f) have multilevel block Toeplitz structure and to the case where 
the test functions F have not bounded support (see [20] and references therein). 

If f is not real-valued, then T n (f) is not Hermitian in general: consequently, the distribution of eigenvalues 
is more involved and (5) cannot be extended in the natural way (see [29]). A very elegant geometric based 
result is due to Tilli [30] and the conclusion is surprisingly simple: a Toeplitz sequence with bounded symbol 
/ will have a canonical eigenvalue distribution in the sense of (1), if the complement of the range of / is 
connected and the range has empty interior. This makes clear that regularity plays no role and this explain 
why this result was not found for many years: researchers were in the wrong direction looking at regularity 
assumptions on the symbol. The same misunderstanding occurred, in minor proportions, for the conditioning 
of a Toeplitz matrix generated by a weakly sectorial symbol [3]: again it is a geometric phenomenon that 
describes the asymptotic behavior of the conditioning and not a regularity property of the symbol. Take 
f(s) = (2 — 2 cos(s)) 10 . Then the minimal eigenvalues of the single-level T n (/) tends to 0 (the infimum of /) 
monotonically and with asymptotic speed dictated by n -20 (notice that 20 is the order of the unique zero 
of /). Exactly the same behavior is proven (with a different constant [15, 3]) if /(s) = (2 — 2cos(s)) 10 h(s) 
where h(s) is any real-valued L°° function with positive infimum: indeed the result is a consequence of how 
the essential range of the nonnegative symbol / “touches” 0 from above and the fact that f(s) is infinitely 
differentiable, as in the case of h(s) = 1, or is discontinuous almost everywhere (a.e.) does play any role. 

2.3 GLT matrix sequences 

For the subsequent analysis, it is convenient to introduce the class of Generalized Locally Toeplitz (GLT) 
sequences that represents at the same time a generalization of Toeplitz sequences and of matrix sequences 
approximating variable coefficient (differential) operators [21]. More in detail, the class of GLT sequences 
can be essentially viewed as a topological closure, both in the matrix side and in the “symbol” side, of linear 
combinations of products of Toeplitz sequences and diagonal sampling matrix sequences: a sampling matrix 
(of level 1) D n (a) of size n is the diagonal matrix containing as j’-th diagonal element a(jh), h a mesh 
parameter, a smooth enough. For our purposes we do not need to introduce the (quite long and involved) 
definition of this class for which we refer to [21, 24]. We just recall the main properties especially those 
which are of interest for our problem. 

A. Any GLT sequence { A n } is uniquely associated to a measurable symbol k(x,s), x € fi Peano-Jordan 

measurable set of R d (space domain), s € Q d (Fourier domain), D = fi x Q d : we write {A n } ~glt k 
and we have {A n } (k,D) and {A n } ~ x (k,D) if A n Hermitian at least for n large enough. 

B. Every Toeplitz sequence generated by f(s) in the sense of (3) is a GLT sequence with k(x,s ) = f(s) 

(Szego-Tyrtyshnikov theory). 

C. Every sequence which is distributed as the zero function in the sense of (1) for the singular value is a 

GLT sequence with k(x, s) = 0. 

D. Every Finite Difference (FD) and Finite Element (FE) equi-spaced approximations of constant coefficient 

PDEs on square regions (any boundary condition) is a GLT sequence with k(x, s) = p(s) for some 
trigonometric polynomial p (Fourier Analysis); 
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E. Any FD, FD discretization of general variable coefficient (system of) PDEs over Q is a GLT sequence. 

In that case k(x,s) is easily identified (Generalized Fourier Analysis): k(x,s) is the principal symbol 
with obvious changes of the Kohn-Nirenberg and Hormander theory for Pseudo-Differential operators. 

The GLT sequences form a *-algebra. More precisely, the GLT sequences are stable under linear com¬ 
binations, product, pseudo-inversion, and adjoint. In fact, if {A„} ~glt k-a and {£?„} ~glt kb, then we 
observe stability under 

F. linear combinations i.e. {aA n + /3B n } ~glt otn A + P^b] 

G. product i.e. {A n B n } ~ G lt k a k b ; 

H. (pseudo)-inversion i.e. {Af} ~glt provided that {A n } is invertible (invertible elements are those 

such that the symbol vanishes at most on a set of zero Lebesgue measure=sparsely vanishing). 

I. adjoint i.e. {A„} ~glt ka is equivalent to {A*} ~glt k* a . 

In the following, we will use essentially properties A., B., C., and the structure of algebra of the GLT 
class. 


3 Identification of the multiplier 

The section is divided in three parts. In the first we discuss in detail the solution to our problem in one 
dimension: as already mentioned, it turns out that the boundedness of cj> is not necessary and only the 
Lebesgue integrability of a related symbol </> is crucial. Subsections 3.2 and 3.3 are devoted to sketch the 
solution for multivariate and matrix-valued multipliers in the case of separable weight functions. Instead 
of giving all the details, we will emphasize what is new in the derivation and the surprise is that the 
multivariate matrix-valued problem does not pose essentially more difficulties than the scalar case in one 
dimension (except, may be, for the notations). 

3.1 The problem in 1 dimension 

Let 0 be a bounded function defined on [—1,1] and let us consider the multiplication operator M[<j>\ : 
L^,([—1,1]) —■> L^,([—1,1]) defined as M[<f>](f) = </>/, w suitable weight function with the usual assumptions. 
It is known that the spectrum is given by the essential range of p (see e.g. [13]): now suppose that we have 
only a finite number of coefficients 

(6) M n [<j>] = ((M[(j)\ej, ei))"Tj 0 , 

with {e^} denoting an orthonormal basis of L 2 W . The question concerns the reconstruction of the multiplier 
<j> from the spectrum of M n [(j)\ in the sense discussed in Subsection 2.1. 

Consider first the case of the Clrebysliev weight w{x) = (1 — a: 2 ) -1 / 2 of first kind and of its basis 
ej(x) = cos(j arccos(a:)). Then 

{M n [<j)]) i :j = (M[(j>]ej,ei) = j <j>(x)ej{x)ei(x) dx = J p(cos(s))cos(js)cos(is)ds, Q = (-7r,7r). 

As a consequence, the matrix M n [f] can be expressed in terms of the Fourier coefficients fk of the function 
/(s) = 7r</>(cos(s)) in the sense of (3) and then (M n [(j)]) i ■ = +fj-i +fi+j +f-i-j)- Taking into account 

that /(s) is even we directly see that fk = f-k for ever k £ Z and therefore 

(Afra^Djj = f\i—j\ + f\i+j\i 

i.e. 

(7) M n [d>\ = T n (f)+H n (f). 

Here the matrix H n (f) = (/|i+j|)".-^ 0 is of Hankel type since its entries are constant along the anti-diagonals. 
Moreover, from [8] we know that the Hankel sequence {H n (f)} is distributed as the zero function over Q 
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in the sense of (1): {H n (f)} is indeed a GLT sequence [21] with symbol equal to zero (see item C.) and 
therefore, since it is bounded in spectral norm, both relationships in (1) hold with g = 0 (the singular value 
part is contained in item A. and C. and the eigenvalue part follows as in Theorem 1.2 of [22]). Consequently, 
the singular value distribution of {M n [(j>] } is decided by the Toeplitz part {T n (f)} and, if <j> is real-valued, the 
same is true for the eigenvalue distribution too: this can be seen directly by using Tyrtyshnikov perturbation 
arguments [31] or, from a more abstract viewpoint, because the GLT class is an algebra that is by A. and 
F. (see also [27]). Therefore the symbol of { M n [<j) ]} = {T„(/)} + {H n (f)} is equal to the one of {T n (f)} i.e. 
/ plus that of {H n (f)} which is zero. 

As a consequence, if the multiplier (f> is real-valued, then we can reconstruct, approximately, <f> from the 
eigenvalues of M n [</>]. In the general case, the desired result depends on the geometric structure of the range 
of (f) and on the Hankel correction: the eigenvalues can be dramatically sensitive even to 1 rank corrections 
(see e.g. [34]). This pathological behavior of the eigenvalues has also good side effects because the effective 
procedure that can be designed (see Subsection 3.1.2) depends exactly on the existence of close sequences 
whose spectral behavior is substantially more regular than Toeplitz sequences (see Remark 3.2). 

The case of the Chebyshev weight of second kind is also very simple to handle thanks to the explicit 
expression of its orthogonal basis elements after the usual change of variable x = cos(s). Indeed we have 
w(x) = (1 — a: 2 ) 1 / 2 and ej{x) = sin((j + 1) arccos(a:))/sin(arccos(x)). Then 

{M n [4>\) i :j = (M[(f>\ej,ei ) = J (f>(x)ej{x)el{x) dx = J 0(cos(s)) sin((j + l)s) sin((z + l)s) ds, Q = (—7r, tt). 

From the latter we infer (. M n \(f>]) i ■ = f\i~j\ + f\i + j + 2 \ with /*, Fourier coefficients of f(s) = 7r^(cos(s)) and 
then M n [(j)] = T n (f ) + H n (f ) with H n {f) being the principal sub-matrix (of size n) made by the last n 
rows and columns of H n+ \(f) and H n {f) as in (7). Therefore a simple interlace argument (see e.g. [1]) 
shows that the corresponding Hankel sequence is distributed as the zero function over Q and then (see [21]) 
{M n [cl)\} = {T„(/)} -I- {#„(/)} has the same GLT symbol as {T n (f) } i.e. / and the conclusion is as before. 

In fact, the above analysis can be generalized using purely linear algebra tools but the result itself is 
known already thanks to Szego (see [26]). For every choice of the weight function the symbol of { M n [f] } is 
always f(s) = 7r^(cos(s)) which is independent of the weight function w. In other words, the finite sections 
of M[(p] with orthogonal polynomials always give more attention to the endpoints of the original interval —1 
and 1 and less attention on the central part of the domain. That behavior is also important for the success 
of many associated numerical methods such as Gaussian quadrature formulae and interpolations schemes at 
the zeros of orthogonal polynomials. 

However, let us give a short look to a sketch of a linear algebra derivation. 

Proposition 3.1 Consider a general weight w with the usual restrictions (nonnegative, with support co¬ 
inciding with [—1,1], with finite Lebesgue integral). Let ej be the j-th orthogonal polynomial. Then the 
following facts hold: 

1. ej(x) = J^i -0 a i c i( x )> a j 7^ °i i~th Chebyshev polynomial of first kind; 

2. E n _i(x) = L n F n -i(x), L n lower triangular invertible matrix, E n -i(x) n-dimensional vector whose 

i-th position, i = 0, 1, is the given by efix) and F n _i(x) n-dimensional vector whose i-th 

position, i = 0,... ,n — 1, is the given by cfix); 

3. M n [(j>] = f\4>(x)w(x)E n _ 1 (x)E0_ 1 (x)dx; 

4■ M n [<j)\ = L n ■ f\ <j>(x)w(x)F n _i (x)F^_ 1 (x) dx ■ ; 

5. M n \j)\ = L n ■ M n [(j)\ ■ L(f, with 4>{x) = 4>(x)w{x)\/l — x 2 and M n [fi] being the n-th finite section of 
M[(j>\ in the case of the Chebyshev weight of first kind; 

6. {M n [4>}} ~<T (f,Q), f(s) = 7T<^(cos(s)) = 7T(/>(cos(s))ie(cos(s))sin(s), Q = (-7r,7r); 

7. {L n } is a GLT sequence with symbol \Jg(s) and {L(f L n } is a GLT sequence with symbol g(s) = 
l/['u>(cos(s)) sin(s)]; 
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8. {L^L n } ~ AjCT (g,Q), g(s) = l/[w(cos(s)) sin(s)], Q = (—tt, tt); 

9. {M n [(j)\} is a GLT sequence with weight f(s) = 7T^(cos(s)) 

10. {M n [(j>]} (/, Q), f(s) = 7r0(cos(s)). 

Discussion on the proof. Item 1. is obvious since every e 3 has degree j and the second item is again 
obvious since e 3 has exactly degree j. Item 3. is a compact rewriting, directly in matrix form, of (6). Item 
4. follows from Item 2. and Item 3. taking into account the linearity of the integral and that L n does not 
depend on x. We now recall that F n _i(x) contains the Chebyshev basis of first kind: therefore, in order to 
interpret the scalar product as the one induced by the Chebyshev weight of first kind, the related multiplier 
has to be seen as <j>(x) = 4>{x)w{x)\/l — x 2 and Item 5. is proved. Further, after the usual change of variable 
x = cos(s), the matrix M n [(j)\ can be written as T n (f) plus H n (f). Moreover, w £ L 1 [—1,1] and therefore 
/ £ L 1 (Q), Q = (—7r,7r). Finally by Theorem 2.2 (which holds for L 1 functions) and by [8], we know that 
the Toeplitz part is distributed as / and the Hankel part as zero, respectively (also Item A., Item B., and 
Item C.). This is enough by Item F. (see Theorem 4.5 and Section 5 in [21] for more details) for deducing 
that {M n [(j )\} ~ (T (/, Q) and Item 6. is proven. The first part of Item 7. is just a technical calculation (see 
[23]) and the second part of the same item and Item 8. follow from the first part and from Item G. and 
Item I. (see also Theorem 4.5 in [21]). Finally, by Item 5. and Item 7., {M n [<j>\ } is a product of two GLT 
sequences, {M n [(j)\} and {L^L n } with symbols 7r</>(cos(s))u;(cos(s)) sin(s) and l/[w(cos(s)) sin(s)], respec¬ 
tively. Therefore, due to the structure of algebra of GLT sequences (again Theorem 5.8 in [21] i.e. Item G.), 
{M n [(j)\} is a GLT sequence with symbol 7r^(cos(s)) = ir<j)(cos(s))w(cos(s)) sin(s)/[w(cos(s)) sin(s)] (Item 
9.) and, finally, {M n [cj)]} ( f,Q ), /(s) = 7r^(cos(s)) that is Item 10., again by Theorem 4.5 in [21] i.e. 

Item A. • 

It is clear that the above proof is really compressed and that some relevant details have to be found in 
other papers (in particular for Item 7.). The point was to show that from purely linear algebra reason¬ 
ings it is possible to treat this kind of problems and sometimes obtaining in a simpler way more general 
information: see e.g. [12] where the analysis of the zero distribution of orthogonal polynomials with vary¬ 
ing coefficients is made by employing GLT arguments, without any regularity assumption except for the 
Lebesgue measurability. We emphasize in addition that the algorithm in the next subsection depends only 
on Items 4., 5., and 6. for which we have given a detailed proof and that item Item 6. is indeed valid 
as long as f £ L 1 {Q). We observe that the latter means that the assumption on the boundedness of the 
multiplier </> is not necessary and can be dropped. More specifically, we can allow <fi to be just Lebesgue 
integrable if we have w(cos(s)) sin(s) £ L°°(Q ): we already encountered examples in this direction and 
namely the Chebyshev weight of first kind for which tu(cos(s)) sin(s) = 1 and that of second kind for which 
w(cos(s)) sin(s) = sin 2 (s). It is clear that there exists a large class of weights of this type. 

3.1.1 Circulant approximation 

We start by describing the circulant class with special attention to its approximation properties with respect 
to Toeplitz matrix sequences. The algebra of circulant matrices is a subclass of Toeplitz matrices to which 
it is not possible to attribute a symbol in the sense of (3) with exception for the identity and for the null 
matrix. In the one-level case (the one discussed so far in this section), they share the algebraic property that 
every row is the forward circular one-step shift of the previous row and where also the notion of “previous” 
has to be intended in a circular way: more precisely, the first row can be seen as the forward circular one-step 
shift of the last row as it is clear from equation (8). The latter nice algebraic feature translates in many 
properties related to circular convolutions. Here we only point out another important characterization in a 
spectral sense. Every circulant matrix of size n can diagonalized by the (unitary) discrete Fourier matrix. 
This means that A n is circulant if and only if A n = F n DF* where I? is a complex diagonal matrix, 

F„ = , k,j = 0,... ,n - 1, 

is the Fourier matrix of size n and X* denotes the complex transpose of X. Moreover, the diagonal matrix 
D has j-th entry given by p n (a^" ) ) with x^ n) = 2nj/n, j = 0,..., n - 1, p n (z) = J2kZ o a kZ k , a 0 ,..., a n _i 
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being the entry of the first column c[l] of A n =circ(a) i.e. 


( 8 ) 


a0 ffln-l 

<22 

<2l 

«1 Oq (Jn-1 


<22 

a 2 





<2n —1 

•n—1 ' ' ' «2 

ai 

<2 0 


Notice that the above eigenvalue formula has also an important computational counterpart since the 
vector d containing the diagonal entries D is equal to F*c[l\ and F* = PF n , with P flip-type permutation 
matrix. As a consequence, the spectral decomposition of any circulant matrix can be recovered in 0(?rlog(n)) 
complex operations via the celebrated FFT (see [33]). We now recall some connections between circulants 
and (one-level) Toeplitz matrix sequences associated to a symbol. 


Definition 3.1 Let C n be the algebra of circulant matrices and let T n (f) be a single level Toeplitz matrix 
associated to the symbol f. Then the following definitions hold. 

• The Strang preconditioner N n (f) associated to T n (f) is the circulant matrix obtained from T n (f) by 
copying the first [n/2] central diagonals with [a;] denoting the rounding of x. In other words, the j-th 
entry of first column c[l] of N n (f), j = 0,..., [n/2] — 1, is exactly the j-th Fourier coefficient aj of f. 

• The optimal preconditioner C n (f) = Opt(T„(/)) is the unique solution of the minimization problem 

(9) min||A-X|| F , A = T n (f ), 

with || • ||f denoting the Frobenius norm i.e. the Euclidean norm of the singular value vector (Schatten 
p norm with p = 2) or, equivalently, the Euclidean norm of n 2 -sized vector obtained by putting in a 
unique vector all the columns of the argument. 


Some remarks are in order. The existence and uniqueness of the Strang or natural preconditioner are 
implicit in the definition itself which clearly indicates an explicit cost-free expression. The existence and 
uniqueness of the optimal preconditioner (see e.g. [5]) follow from the strict convexity of the Frobenius norm 
that implies the existence and uniqueness of the minimizer from a given convex closed set. We are in a finite 
dimensional setting and clearly the linear space of the circulants C n is closed and convex. 

Finally, the optimal approximation admits an easy to derive and very interesting representation since 

(10) Opt(A) = F n diag (F* AF n ) F*, 

where A ia a generic complex square matrix and the operator diag applied to any square matrix X gives the 
diagonal matrix whose diagonal entries coincide with those of X. Moreover, if A = T n (f ) then 

(11) Opt(A) = C n (f) = circ(a), a* = ^ + ifi- n ) (i = 0 ,...,\n - lj). 

In the nest proposition we discuss the spectral properties of these matrix approximations by focusing on 
the relationships with the related approximation of the symbol. 

Proposition 3.2 Let f € L 1 (Q), Q = (—n, it), and let us consider N n (f) and C n (f) be the Strang and 
optimal approximations ofT n (f), respectively. Then the following facts hold: 

1. The Strang preconditioner N n (f) has eigenvalues lF n >[f}(x^), j = 0,... ,n — 1, where n' = [n/2] — 1, 
and P q [f] is the Fourier sum of degree q of f (see [5]j. 

2. In the general case where f € L X {Q) and it is not smooth, anything can happen: N n (f) definitely 

singidar or indefinite even ifT n (f ) is positive definite for every n, N n (f) collectively unbounded even 
if \\T n (f)\\ < H/lloo for every n, {N n (f)} clustered at infinity even if {T n (f)} (/, Q). 

3. If f belongs to the Dini-Lipschitz class and is 2ir-periodic, then the eigenvalues of N n (f) will reconstruct 
f in uniform norm. 
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f. The optimal preconditioner C n (f ) = Opt(T„(/)) has eigenvalues C n -i[f](Xj), j = 0,..., n— 1, where 
C q [f} = is the Cesaro sum of degree q of f (see [ 17 }). 

5. If f is continuous and 2ir-periodic, then the eigenvalues of C n (f) will reconstruct f in uniform norm. 

6. ||Opt(7l)||* < ||7l||* for every unitarily invariant norm and in particular ||C' n (/)|| p < ||T n (/)||j, for 
every Schatten p norm, p> 1 (see Theorem 2.1, item 6, in [7]). 

7. Iff is L°°(Q), then \\C n (f)\\ < ||/||oo and, if f £ L P (Q) then \\C n (f)\\ p < £ f Q \f(s)\ p ds. 

8. {C n (f)} distributes as ( f,Q ) both in the sense of the eigenvalues and singular values. 

9. With the notation of Proposition 3.1, {Opt(M„[^>])} disti'ibutes as ( f,Q ) both in the sense of the 
eigenvalues and singidar values with f(s) = 7r</>(cos(s))w(cos(s)) sin(s). 

Proof. Items 1., 4., 6. can be found in the relevant literature. Item 3. is a direct consequence of 
the fact that the Lebesgue constant of the Fourier sum is asymptotic (up to a multiplicative constant) to 
log(?r) and therefore the Fourier sum has to converge to / since the modulus of continuity of / satisfies 
uif{l/n) = o( 1/ log(n)) for every / in the Dini-Lipschitz class. Item 2. is a nice application of known facts. 
The example of Du Bois-Raymond is a nonnegative function / £ L°°(Q) with unbounded, highly oscillating 
Fourier sum (see e.g. [2]). Clearly the matrix N n (f) is unbounded and definitely indefinite while T n (f) is 
positive definite and uniformly bounded in spectral norm by ||/||oo (for the Toeplitz part see e.g. [25] where 
also the tools for proving item 6 of Theorem 2.1 in [7] can be found). For finding an example where {iV„(/)} 
clustered at infinity even if {T n (/)} (f,Q), it is enough to use the example of Kolmogorov (see e.g. [2]): 

the function belongs to L l (Q), but it is not in L 2 (Q) and has a Fourier sum diverging everywhere so that the 
eigenvalues of N n {f) collectively explode, but thanks to Theorem 2.2 it is still true that { T n (f)} ~ (7 (/, Q). 
Item 5. is trivial since (thanks e.g. to the beautiful theory by Korovkin) it is well known that the Cesaro 
sum of any continuous function converges uniformly to /. By [25] we know that ||T„(/)|| < ||/||oo whenever 
/ € L°°(Q) and ||T„(/)||p < ^ Jq \.f(s)\ p ds whenever / € L P (Q) with p > 1: as a consequence, Item 7. 
follows from Item 6. 

Concerning Item 8. we remark it has been proved that for every f £ L 1 (Q), {T n (/)} (f,Q) and 

{T n (f)} ~a (/, Q) if / is real-valued (see [18]). We then need only to prove that the distribution results 
stands for the eigenvalues as well even for complex-valued symbols (notice that the latter is not trivial since 
it does not hold in general in the Toeplitz case as observed by Morrison in [13]). 

We want to prove that 

lim S x (F,C n (f)) = ^~ [ F(f(s))ds 

for every / £ L 1 (Q), for every F continuous with bounded support in C. First we observe that the claim can 
be reduced to the case of F Lipschitz continuous with bounded support in C. In fact for every G continuous 
with bounded support in C, for every e > 0, we can find G e Lipschitz continuous with bounded support 
such that | G(z) — G e (z)\ < e for every z £ C (notice that in general we cannot take G e polynomial due to 
the obstruction given by the Mergelyan theorem (for a proof see [14])). 

Now by Item 5. the claim is already proven if / is continuous and 27r-periodic (notice that in the Toeplitz 
case this is again false in general with elementary polynomial examples). Therefore for every / £ L 1 (Q), for 
every e > 0, we consider f e continuous and 27r-periodic such that ||/ — /e||i,i(Q) < 27re so that 

Y [ nf(s))ds~Y f F(f e (s))ds f \F(f(s)) — F(f e (s))\ ds < Me 

J Q J Q J Q 

with M being the Lipschitz constant of F. Moreover, by the same argument we have, 

-t n 1 n 

|S A (F,C„(/))-E A (F,C'„(/ e ))| < - Y^in^Cnim-FiXjiCMem < 

" 3 = 1 " 3 = 1 

and, since the circulants form an algebra and the operator C n (-) is linear, we have 

1 n 

|E x(F,C n (f)) - E x(F,C n (fe))\ < M- £ |A 0 {C n {f - / e ))|. 

n ■“ 
j=i 
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But the singular values of any circulant matrix are exactly the moduli of its eigenvalues since every circulant 
is also normal. Therefore, by Item 7., we have 

|£ A (F,C„(/))-£ A (F,C„(/ e ))| < ^\\C n (f-fe)\\i<Me 
and the proof is concluded since e is arbitrary. 

We conclude with the proof of Item 9. By Item 5. and Item 6. of Proposition 3.1, we have 
M n [(j>] = T n (f) + H n (f) and f(s) = 7T^(cos(s))w(cos(s)) sin(s). Therefore by linearity of the operator Opt(-) 
we deduce Opt (M„ [(/>]) = Opt(T„(/)) + Opt = C n (f) + Opt (H n (f)). Now, by the previous item, 
{C n (fJ} distributes as / over Q both in the sense of the eigenvalues and singular values. Moreover, by [8], 
||R„(/)||i = o(n ) and therefore by Item 6. ||Opt(i/ n (/))||i = o(n). Furthermore, from Lemma 2.1 we deduce 
{Opt (H n (f))j (0, Q), {Opt(M n [0])} (f,Q) and {Opt(#„(/))} ~ A (0,Q), {Opt(M„[0]){ ~ A (f,Q) 

if <p is real-valued. Finally, for the complex-valued case when considering the distribution in the eigenvalue 
sense, the proof is as in the preceding item. • 


Remark 3.1 Item 1. in the above proposition has an interesting consequence. Take / £ L°°(Q) and 
consider N n (f). Since the entries of N n (f) contain exactly the same coefficients as T n r(f) with every Fourier 
coefficient counted 2 n' times it follows that ||./V T! (/)||| = ^ r || F „/[/]||| 2 < < «||/||^o- Therefore, by 

the spectral decomposition of N n (f) in Item 1., it follows 

n —1 

ll^(/)ll2 = El^'[/K^ ri) )l 2 <n||/|lL- 

3=0 

Consequently, the cardinality of the set of indices j such that T n ‘ [/] (xj ) is unbounded as n tends to infinity 

has to be o(n) and the infinity norm of F n '[f] over the grid-sequence {Xj} n is at most 0(y/n). This means 
that the set of grid points in which the Fourier sum can diverge is negligible and more precisely its cardinality 
is o[n). Taking into account the possible maximal growth of a polynomial of degree n' = [n/ 2] — 1, it follows 
that the set where the Fourier sum can diverge in [— tt, tt] has to be of zero Lebesgue measure and this is a 
linear algebra version of a Car lesson-type result (see e.g. [2]). 

Remark 3.2 In [13], the author observed that Toeplitz sequences are unable to reconstruct /, in general, 
if / is complex-valued. Tilli gave a precise answer by characterizing the cases where this reconstruction is 
just impossible. In Item 8., we proved that a special circulant approximation of T„(/) is indeed able to 
reconstruct the symbol / in the maximal generality that is for / £ L 1 (Q): moreover, by Item 5., if / is 
also continuous and 27r-periodic, the reconstruction can be performed in a strong sense i.e. in uniform norm. 
This is confirmation of the great stability of the considered approximation which has two reasons: the first 
is the normality of circulants (in contrast with Toeplitz matrices generated by complex-valued symbol which 
can be of maximal non-normality as any Jordan block), the second is the stability of the Frobenius optimal 
approximation which has to be related to the stability of Linear Positive Operators. What we will discuss 
in the next subsection is interesting, because it shows that, under mild assumptions, the problem of the 
identification and reconstruction of the multiplier (/> can be reduced also to the Frobenius optimal circulant 
approximation of a Toeplitz matrix generated by a L l (Q) symbol: the theoretical basis relies on Proposition 

3.2 and especially Items 5., 8., and 9. 

3.1.2 Circulant based algorithms 

We only suppose to know the coefficients of M n \cj>] and the weight w with the related matrix L n and unknown 
4> (the case where (f> is known with unknown weight w leads to a different problem). The algorithm is heavily 
related to the analysis in Proposition 3.1 and in Proposition 3.2. It can be roughly sketched as follows: 

1. Form M n [(j)\ and from L n (known when the weight w is known), compute X n = M n [(j>] with <j>(x) = 
4>{x)w{x)\/l — x 2 ; 

2. compute C n the Frobenius optimal approximation of X n ', 
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3. compute the eigenvalues of C n by FFT (storing also the index of the related eigenvectors); 

4. reconstruct the function f(s) = n4>(cos(s))w(cos(s)) sin(s) and therefore dividing by 7ru;(cos(s)) sin(s) 
reconstruct /(s) = ^(cos(s)), s G Q i.e. <j>(x), x G [—1,1]. 

Indeed the correctness of the above procedure is based on the last item of Proposition 3.2, since X n = 
M n [<f)]) and C n = Opt(M n [0]) (see also Item 5. and Item 8.). 

We observe that the matrix X n in view of (7) contains an Hankel part which represents a disturbance. 
Therefore, also in order to exploit the computationally convenient formula (11), we can eliminate this part. 
The argument is a trivial application of the Riemann-Lebesgue Lemma (see [14]): indeed, instead of X n = 
M n [<j>\ = T n (f) + H n (f) we would like to consider the matrix T n (f) only. Unfortunately, the matrix T n (f) is 
unknown (only the entries of the whole matrix X n are available) and we will approximate it by the Toeplitz 
matrix T n constructed according to the following idea. We have (X n ) n<n = {T n (f)) n<n + ( H n (f)) n<n = 
fo + f 2 n ~ fo since, by the Riemann-Lebesgue Lemma, / 2n is infinitesimal: we set (T n )j,j = fo + f 2 n , 
j = 1, . . • , U. We observe that (W n ) n _i >n + (^ln)n,n — 1 = (^n(/))n-l,n A {T n (f T (-^n(/))n— l,n + 

(7L„(/)) n , n _ 1 = /—i + fi + 2/ 2 n-i- Since / is even we have f-j = fj for all j G Z and therefore ((X n ) n _i >n + 
{X n ) n<n -i)/2 = fi + / 2 n—l ~ fi since, by the Riemann-Lebesgue Lemma, / 2n -i is infinitesimal: we set 
1 = {T n )j-i,jfi + hn-x, j = 2,...,n. We proceed by considering (X„) n _ 2j „ + (X n ) ra _i j „_i + 
{X n ) n<n -2 = f —2 + fo + f 2 + 3/ 2n _i. Now we already computed an approximation of /o and we know that 
/_ 2 = / 2 . Therefore we can compute / 2 within an infinitesimal error since 

j-2 = {T n )j-2,j = [((X n ) n _2,n + (^n)n-l,ra-l + (- 5 fn)n,ra-2) — (^n)n,ra] /2 = /2 + (3/ 2n -l — /2n)/2 ~ / 2 , 

j = 3,... ,n. The procedure can be continued in a similar way by obtaining every entry of T„(/) i.e. every 
Fourier coefficient fj, |j| < n — 1, within an infinitesimal approximation error. 

The new algorithm can be written as follows where the third step is obtained by the previous reasoning. 

1. Form M n [(f>] and from L n (known when the weight w is known), compute X n = M n [(j>\ with 4>{x) = 
4>(x)w(x)\/l — x 2 ; 

2. compute T n approximation of T„(/); 

3. compute C n the Frobenius optimal approximation of T„; 

4. compute the eigenvalues of C n by FFT (storing also the index of the related eigenvectors); 

5. reconstruct the function f(s) = n(/)(cos(s))w(cos(s)) sin(s) and therefore dividing by 7rw(cos(s)) sin(s) 
reconstruct /(s) = <(>(cos(s)), s G Q i.e. 4>(x ), x G [—1,1]. 

The correctness of the algorithm depends entirely on the fact that {T„} distributes as (/, Q) in the sense 
of the singular values. We know {T n (f)} ( f,Q ) and | (T n — T n {f))j k \ tends to zero for every ( j,k ) 

(more precisely we have \(T n — T n (f))j^ k \ = 0(/„)). Unfortunately, the second relation does not imply that 
||T„ — T n (f)\\ p = o(n 1 / p ) for some p G [1, oo] and therefore we are not allowed to use Lemma 2.1. In fact, let 
us consider the following example. Assume that 

T n - T n (f) = e n V^F n , e n > 0 infinitesimal. 

Then every entry has modulus e„ but every eigenvalue has modulus equal to \fnt n and therefore {T„ — T n (f)} 
distributes as (0,Q) only if e n = o(?r -1 / 2 ). Indeed, defining e n the maximal values of |(T„ — T n (f))j , k |, we 
have 

n —1 n —1 

l|r„-T n (/)|||= 5] \(T n — T n (f))j tk \ 2 < J2 4 = 4« 2 - 

j,k —0 j,k=0 

Therefore, for p = 2 we obtain ||T„ — T n (f)\\ p = o(n 1 / 2 ) if /„ = o(n -1 / 2 ). As a consequence, in order to use 
the second algorithm, we should have more information on the symbol /: for instance if / is 27r-periodic and 
fc-times continuously differentiable, k > 1, then /„ = o(n~ k ) and therefore e n = o{n~ k ) so that we can use 
safely the second algorithm and, in addition, if k is large then T n is a very good approximation of T„(/). In 


12 



conclusion, for a smooth symbol i.e. for large fc, the reconstruction provided by the latter algorithm could 
be better than the one given by the first. Of course, one should have this information a priori or, possibly, 
one should use the first algorithm for obtaining a guess: if the result looks like a smooth function (under the 
assumptions of Item 5. in Proposition 3.2, we recall that the approximation is in uniform norm), then the 
second algorithm could be employed for improving the quality of the reconstruction. 

3.2 Generalization: the multivariate separable case 

Let 0 be a bounded function defined on [—1,1]“ and let us consider the multiplication operator M [d] : 
L 2 ([-1, l] d ) —► L 2 W {{- l,l] d ) defined as M[<j>](f) = (j>f with w(x) = Wi(xi)w 2 (x 2 ) ■ ■ ■ Wd(xd), d-variate 
weight with Wj standard univariate weight function. As in the single-variate case the spectrum coincides 
with the essential range of <f>. Consider to have a finite number of coefficients 

(12) M„M=((MMe i ,e l ))”: = 1 0 T , 

n = (ni,..., rid), * = (ii, ■ • •, id), j = (j l, • • • ,jd), 1 vector of all ones as in Subsection 2.2, with {e^} denoting 
an orthonormal basis of L^([— 1, l] d ) defined by ej(x) = (aq) • • • Zj d (xd) with { ej k } orthonormal basis of 
Lw k ([— 1, l])i k = 1,... ,d. The question is again the reconstruction of the multiplier <§> from the spectra of 
M n [<j>] in the sense discussed in Subsection 2.1. 

Take the case of the d-level Chebysliev weight of first kind w(x) = nti(i — x 2 k ) 1 / 2 and of its basis 
e j(x ) = ej k (®i) • • • ej d (xd), ej k ( Xk ) = cos (jk arccos(xfc)), k = 1,..., d. Then with the usual change of variable 
Xk = cos(sfc), k = 1,..., d, we have 


P d 

(M n [4>])i j = (M[(j>]ej,ei) = / ^(cos(si),..., cos(s d )) TT cos(j fc s fc ) cos(z fc s fc ) ds, Q = (—K,n). 

J Q d kJi 

As a consequence, the matrix M n \(j>\ can be expressed in terms of the d -indexed Fourier coefficients fj of the 
function f(s) = (7r) d ())(cos(si),..., cos (sd)) and then, in d-index notation and taking into account that f(s) 
is even with respect to every variable Sj, we find 

{M n [<t>])ij = f\i—j\ + f\i+j\ 

and therefore 

(13) M n [<f,]=T n {f)+H n {f). 

Once we arrive here, the rest is a straightforward generalization of the univariate case since the result on 
Hankel matrices (see [8]) are directly stated in an arbitrary number of dimensions. Theorem 2.2 is in d 
dimensions and the same applies to the results on the GLT class whose definition is inherently d-dimensional 
(see Subsection 2.3 and [21]). Furthermore, formulae (10)—(11) stand unchanged (d-indices in place on simple 
indices, F n = F ni (g) • • • (g> F nd , n in the denominator of (11) replaced by N(n)) and Proposition 3.2 is again 
unchanged. Therefore also the algorithms can be described verbatim and therefore the optimal circulant 
approximation of M n [(j>\ ss T n (f ), T n « T n (f) will reconstruct with infinitesimal error the Cesaro sum of the 
function / and therefore of (f>. 

3.3 Generalization: the multivariate separable matrix-valued case 

Let (j) be a bounded function defined on [—1, l] d and having values in the space C pxq and let us consider 
the multiplication operator : A 2 ([—1, l] d , C qxr ) —> L 2 ,([— 1, l] d , C pxr ) defined as M[(j>]{f) = (j>f being 

C pxr with / being C qxr and with w(x) = w±(xi)w 2 (x 2 ) ■ ■ ■ Wd{xd) as in the previous subsection. In the 
present matrix-valued setting, it is less obvious to refer to the spectrum of the continuous operator and this 
is true in the discrete as well since the resulting sections are not square matrices. However we can give again 
a meaning passing to the “absolute value” of the operator (see [1]) i.e. the square root of the adjoint times 
the operator itself. In the discrete we are talking of the singular values and in the operator case we are 
talking of the singular values of the multiplier (j). In any case, as we will see in the rest of the derivation, we 
will able to reconstruct <f> (or better its Cesaro sum) and therefore its singular values. 
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It should be observed that the present multiplication operator can be written as a vector whose entries 
are sum of scalar multiplication operators. More precisely we have 

/ q \P’ r 
<t>f = ( ] 

Vt=i / s= i jZ= i 

with 4> s ,tft,z defining a scalar multiplication operator on L^([—1, l] d ). Therefore we can represent M[cj)\ as 

p q 

S=1i=l 

where E(s,t) denotes the p x q matrix being 1 at position ( s,t ) and zero otherwise: notice that {E(s,t)} 

forms the canonical basis C pxq . Therefore, if we consider a finite number of coefficients we have 

M n[<t>\ = ({M[(t>\e j ,e i ))™J* 0 

p q \ n — l 

(^ J '^2,M[(t>s,tE{s,t)}e j ,e i ) j 

S=1 t—1 / i^j — Q 

= ii 

S = 1 t -1 

n = (?ii,..., rid), i = (*i, • • •, id), j = (ji, ■ • • ,jd), 1 vector of all ones, with {ey} denoting an orthonormal 
basis of L 2 W {\— 1, l] d ) as in the latter section. In other words is a multilevel block matrix where the 

size of each block is dictated by the multiplier </> and more specifically we can write 

T (M n [<j >••• {M n [<t> ltq ]) itj 

(14) (M n ^]) Sit = {{M[<j) Stt \ej,ei))) s t g=i (=1 = i i 

The above block expression makes clear that the reconstruction of every single entry ^> Sjt can be done 
exactly as in the scalar-valued case. More precisely the next scheme can be followed. 

• The reconstruction of the entry can be done via the same algorithms proposed in Subsection 3.1.2, 

by extracting from M n [cf)] only the entries of according to (14). 

• If one is interested in a global reconstruction e.g. of the singular values of the multiplier </>, then some 
appropriate new tools have to be introduced: this issue is briefly considered in the next section. 

3.3.1 Toeplitz sequences generated by matrix-valued sumbols 

Let / be a d variate p x q matrix-valued (Lebesgue) integrable function, defined over the hypercube Q d , 
with Q = (—7r,7r) and d > 1. Here the Lebesgue integrability means that every entry of the symbol is a 
standard complex-valued L 1 function. With respect to the notion of matrix-valued symbol, we observe that 
the definition of the coefficients in (3) is formally identical: the only obvious difference is that every fj will be 
a matrix of size p x q. However, the formulae in (1) and consequently Theorem 2.2 do not make sense since 
the eigenvalues and singular values are still scalar while the function to be integrated in the right hand-side 
is p x q matrix-valued. On the other hand, a generalization of that results exists and is quite natural. We 
write that {A n } {f,K,p) and {A n } ~ A if 

(15) lim Y, a (F,A n ) = f jtv[F(\f\)]dn, 

n^oc n(K ) J K l 

1 f 1 

lim T,\(F, A n ) = / -tr [F(f)]dp, l = p = q, f Hermitian-valued, 

rwoo p(R ) J K l 
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respectively, with l being the minimum between p and q, |/| = (f H Z) 1 / 2 and tr[g] = ^2j^j{g), ^j(g), 
j = 1, ... ,1, being the eigenvalues of g. With these notations, we have that any Toeplitz sequence {T„(/)} 
with matrix-valued / € L l (Q d ) (which is equivalent to require that maximal singular value of / is in L 1 {Q d )) 
is such that (15) holds with g being the Lebesgue measure, I\ = Q d , and f being the symbol (see [29, 28, 16]). 
Notice that if g = |/| then tr[g] = JT aj(f) and therefore (15) represents a natural generalization of (1). 

The question is about the reconstruction of the multiplier <j> from the finite section M n \</>]. 

Take the case of the d-level Clrebysliev weight of first kind w(x) = nti(i — x 2 j k ) 1 / 2 and of its basis 
e j(x ) = ej 1 {x\) ■ ■ ■ ej d (xd), ej k (xk) = cos {jk arccos(a;fc)), k = 1,..., d. Then with the usual change of variable 
Xk = cos(sfc), k = 1,..., d, we have 


P d 

(M n [<f>]) ij = (M[<)>]e j ,ei)= / ^(cos(si),..., cos(s d )) TT cos(j fc s fc ) cos(z fc s fc ) ds, Q = {— 

•* Q d fc =i 

As a consequence, the matrix M n \</>] can be expressed in terms of the d -indexed Fourier coefficients fj of the 
function f(s) = (7r) d (/)(cos(si),..., cos(sd)) and then, in d-index notation, we find 

{Mn[4>])ij = f\i-j\ + f\i+j\- 

Taking into account that f(s) is even we directly see that 
(16) M n [<f>]=T n (f)+H n {f). 

Once we arrive here, the rest is a generalization of the multivariate case since the result on Hankel matrices 
(see [8]) are directly stated in an arbitrary number of dimensions with blocks of fixed dimension. Theorem 
2.2 is replaced by the block relation (15) and GLT class has a natural block generalization. Furthermore, 
the formula (11) stands unchanged (d-indices and block coefficients in place on simple indices and scalar 
coefficients) and Proposition 3.2 is again unchanged. Therefore also the algorithms can be described verbatim 
and therefore the optimal circulant approximation of M n [(f>] ss T n (f), T n ss T n (f) will reconstruct with 
infinitesimal error the Cesaro sum of the function / and therefore of (j). Furthermore, by using (15), it is 
possible to reconstruct information on the singular values of / and then of <j>. More precisely, given seR, 
s > 0, there is constructive test, analogous to those in e) and f) of Theorem 2.1, that, starting from the 
singular values of the matrix M n [<j>] « T n (/) (or T n ss T n (/)), tells one if s belongs to the union of the 
(essential) ranges of the singular values of / (and therefore of <p): for details and numerical experiments in 
this block setting see [18]. 

4 Conclusions 

In this paper we have considered the reconstruction of the scalar-valued/multivariate/block-valued multi¬ 
pliers of proper multiplication operators through structured linear algebra tools. Further generalizations 
could be considered as a general compact domain or non-separable weight functions. However, we think that 
the results presented in this note clearly show the utility of purely asymptotic linear algebra tools in the 
considered type of problems in approximation theory. 
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